\documentclass[reprint,amsmath,amssymb,aps,prx,superscriptaddress]{revtex4-2}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{hyperref}

% \newcommand{\supp}{See Supplementary Material [url] for detailed proofs for mathematical conclusions in the main text, numerical results of quasi-1D and strip winding numbers, and the derivation of SGBZs for 2D HN model, which includes Ref.\cite{hatcherAlgebraicTopology2001}}

\begin{document}
\title{Supplementary Materials for \\ Non-Bloch Band Theory for 2D Geometry-Dependent Non-Hermitian Skin
Effect}
\author{Chenyang Wang}
\affiliation{State Key Laboratory of Low-Dimensional Quantum Physics, Department
of Physics, Tsinghua University, Beijing 100084, China}
\author{Jinghui Pi}
\affiliation{The Chinese University of Hong Kong Shenzhen Research Institute, 518057
Shenzhen, China}
\author{Qinxin Liu}
\affiliation{State Key Laboratory of Low-Dimensional Quantum Physics, Department
of Physics, Tsinghua University, Beijing 100084, China}
\author{Yaohua Li}
\affiliation{State Key Laboratory of Low-Dimensional Quantum Physics, Department
of Physics, Tsinghua University, Beijing 100084, China}
\author{Yong-Chun Liu}
\email{ycliu@tsinghua.edu.cn}

\affiliation{State Key Laboratory of Low-Dimensional Quantum Physics, Department
of Physics, Tsinghua University, Beijing 100084, China}
\affiliation{Frontier Science Center for Quantum Information, Beijing 100084, China}

\maketitle
\onecolumngrid
\tableofcontents{}



\setcounter{section}{0} 
\global\long\def\thesection{S\arabic{section}}%
 \setcounter{figure}{0} 
\global\long\def\thefigure{S\arabic{figure}}%
 \setcounter{equation}{0} 
\global\long\def\theequation{S\arabic{equation}}%

% \appendix

\section{Density of Quasi-1D Zeros on PMGBZ}

\label{app:DOS-on-1D-GBZ}

In the main text, we have given the relation between the density of
quasi-1D zeros and the relative phases of PMGBZ pairs without proof,
that is, given a reference energy $E$, the number of $\beta_{1}$-zeros
of $F(E,\beta_{1})=0$ located on a segment of PMGBZ is proportional
to the change of relative phases of PMGBZ pairs at both endpoints.
In this section, we will derive the relation with a modified method
of Ref. \cite{yokomizoNonBlochBandTheory2019}. It should be noted
that the reference energy $E$ is viewed as a fixed parameter in this
section.

In general, consider the hybrid Hamiltonian in the following form,
\begin{align}
\mathcal{H}\left(\beta_{1}\right)= & \sum_{\mu,\nu=1}^{m}\sum_{r_{2}=1}^{L_{2}}\mathcal{T}_{0,\mu,\nu}\left(\beta_{1}\right)c_{\beta_{1},r_{2},\mu}^{\dagger}c_{\beta_{1},r_{2},\nu}\nonumber \\
 & +\sum_{t_{2}=1}^{t_{R2}}\sum_{r_{2}=1}^{L_{2}-t_{R2}}\sum_{\mu,\nu=1}^{m}\left[\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)c_{\beta_{1},r_{2}+t_{2},\mu}^{\dagger}c_{\beta_{1},r_{2},\nu}+\mathcal{T}_{-t_{2},\mu,\nu}\left(\beta_{1}\right)c_{\beta_{1},r_{2},\mu}^{\dagger}c_{\beta_{1},r_{2}+t_{2},\nu}\right],\label{seq:quasi-1D-Hamiltonian}
\end{align}
where $\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)$ is the
$\beta_{1}$-dependent coupling coefficient, which is, 
\begin{equation}
\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)=\sum_{t_{1}=-t_{R1}}^{t_{R1}}\mathcal{T}_{t_{1},t_{2},\mu,\nu}\beta_{1}^{-t_{1}}.\label{seq:T-beta1}
\end{equation}
The non-Bloch waves in the minor direction is determined by the 2D
non-Bloch Hamiltonian, 
\begin{equation}
h_{\mu\nu}\left(\beta_{1},\beta_{2}\right)=\sum_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)\beta_{2}^{-t_{2}},\label{seq:2D-non-Bloch-H}
\end{equation}
and the corresponding ChP reads, 
\begin{align}
f\left(E,\beta_{1},\beta_{2}\right) & \equiv\det\left[E-h\left(\beta_{1},\beta_{2}\right)\right],\nonumber \\
 & =\begin{vmatrix}E-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},1,1}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & -\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},1,2}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & \cdots & -\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},1,m}\left(\beta_{1}\right)\beta_{2}^{-t_{2}}\\
-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},2,1}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & E-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},2,2}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & \cdots & -\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},2,m}\left(\beta_{1}\right)\beta_{2}^{-t_{2}}\\
\vdots & \vdots & \ddots & \vdots\\
-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},m,1}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & -\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},m,2}\left(\beta_{1}\right)\beta_{2}^{-t_{2}} & \cdots & E-\sum\limits_{t_{2}=-t_{R2}}^{t_{R2}}\mathcal{T}_{t_{2},m,m}\left(\beta_{1}\right)\beta_{2}^{-t_{2}}
\end{vmatrix}.\label{seq:det-2D-ChP}
\end{align}
According to Eq. (\ref{seq:det-2D-ChP}), when $\mathcal{T}_{t_{2},\mu,\nu}\left(\beta_{1}\right)$
are non-zero for any $\mu,\nu=1,2,\dots,m$, the highest and lowest
degrees of $\beta_{2}$ in $f\left(E,\beta_{1},\beta_{2}\right)$
are $-M_{2}=-mt_{R2}$ and $N_{2}=mt_{R2}$, respectively. Therefore,
for given values of $E$ and $\beta_{1}$, there are $M_{2}+N_{2}=2mt_{R2}$
$\beta_{2}$-zeros for the eigenvalue equation $f\left(E,\beta_{1},\beta_{2}\right)=0$.
We sort the zeros by $\left|\beta_{2}^{(1)}\right|\leq\left|\beta_{2}^{(2)}\right|\leq\cdots\leq\left|\beta_{2}^{(2mt_{R2})}\right|$.
The corresponding non-Bloch waves in the reciprocal space can be calculated
by, 
\begin{equation}
h\left(\beta_{1},\beta_{2}^{(j)}\right)\tilde{\phi}^{(j)}=E\tilde{\phi}^{(j)},\label{seq:non-Bloch-wave-reciprocal}
\end{equation}
where $\tilde{\phi}^{(j)}=\left(\tilde{\phi}_{1}^{(j)},\tilde{\phi}_{2}^{(j)},\dots,\tilde{\phi}_{m}^{(j)}\right)\in\mathbb{C}^{m}$,
and the corresponding real-space expression reads, 
\begin{align}
\phi_{\mu}^{(j)}\left(\beta_{1},r_{2}\right) & \equiv\langle\beta_{1},r_{2},\mu|\phi^{(j)}\rangle=\tilde{\phi}_{\mu}^{(j)}\left(\beta_{2}^{(j)}\right)^{r_{2}},
\end{align}
where $\left|\beta_{1},r_{2},\mu\right>\equiv c_{\beta_{1},r_{2},\mu}^{\dagger}\left|0\right>,r_{2}=1,2,\dots,L_{2},\mu=1,2,\dots,m$
is the single-particle basis within a supercell.

Next, we assume that the OBC eigenstate of $\mathcal{H}\left(\beta_{1}\right)$
has the form of Eq. (\ref{seq:formal-OBC-eigenstate}), 
\begin{equation}
\left|\psi\right>=\sum_{j=1}^{2mt_{R2}}C_{j}\left|\phi^{(j)}\right>.\label{seq:formal-OBC-eigenstate}
\end{equation}
Substituting Eqs. (\ref{seq:non-Bloch-wave-reciprocal}-\ref{seq:formal-OBC-eigenstate})
into the eigenvalue equation $\mathcal{H}\left(\beta_{1}\right)\left|\psi\right>=E\left|\psi\right>$,
the equations of $C_{j}$ read, 
\begin{align}
\psi_{\mu}\left(-l\right) & \equiv\sum_{j=1}^{2mt_{R2}}C_{j}\tilde{\phi}_{\mu}^{(j)}\left(\beta_{2}^{(j)}\right)^{-l}=0,\label{seq:eq-C-1}\\
\psi_{\mu}\left(L_{2}+1+l\right) & \equiv\sum_{j=1}^{2mt_{R2}}C_{j}\tilde{\phi}_{\mu}^{(j)}\left(\beta_{2}^{(j)}\right)^{L_{2}+1+l}=0,\label{seq:eq-C-2}
\end{align}
for $l=0,1,\dots,t_{R2}-1$ and $\mu=0,1,\dots,m$. Equations (\ref{seq:eq-C-1})
and (\ref{seq:eq-C-2}) are homogeneous linear equations of $C_{j},j=1,2,\dots,mt_{R2}$,
so that the condition for non-zero solutions of $\left|\psi\right>$
requires the determinant of the coefficients to vanish, i.e., 
\begin{equation}
\begin{vmatrix}\tilde{\phi}^{(1)} & \tilde{\phi}^{(2)} & \cdots & \tilde{\phi}^{(2mt_{R2})}\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{-1} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{-1} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{-1}\\
\vdots & \vdots &  & \vdots\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{-t_{R2}+1} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{-t_{R2}+1} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{-t_{R2}+1}\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{L_{2}+1} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{L_{2}+1} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{L_{2}+1}\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{L_{2}+2} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{L_{2}+2} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{L_{2}+2}\\
\vdots & \vdots &  & \vdots\\
\tilde{\phi}^{(1)}\left(\beta_{2}^{(1)}\right)^{L_{2}+t_{R2}} & \tilde{\phi}^{(2)}\left(\beta_{2}^{(2)}\right)^{L_{2}+t_{R2}} & \cdots & \tilde{\phi}^{(2mt_{R2})}\left(\beta_{2}^{(2mt_{R2})}\right)^{L_{2}+t_{R2}}
\end{vmatrix}=0.\label{seq:det-of-C}
\end{equation}
It is noted that $\tilde{\phi}^{(j)}$ in Eq. (\ref{seq:det-of-C})
is an $m$-vector, so the determinant is a $2mt_{R2}\times2mt_{R2}$
determinant. According to the definition of $\tilde{\phi}^{(j)}$,
$\tilde{\phi}^{(j)}$ is independent of $L_{2}$, so the determinant
in Eq. (\ref{seq:det-of-C}) can be expanded in the general form of
Eq. (\ref{seq:new-eq-C}), 
\begin{equation}
\sum_{1\leq j_{1}<j_{2}<\cdots<j_{M_{2}}\leq2M_{2}}\left(\beta_{2}^{(j_{1})}\beta_{2}^{(j_{2})}\cdots\beta_{2}^{(j_{M_{2}})}\right)^{L_{2}}g_{j_{1}j_{2}\dots j_{M_{2}}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)=0,\label{seq:new-eq-C}
\end{equation}
where the functions $g_{j_{1}j_{2}\dots j_{M_{2}}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)$
are independent of $L_{2}$. When $L_{2}$ tends to infinity, by dividing
both sides of Eq. (\ref{seq:new-eq-C}) by $\left(\beta_{2}^{(M_{2}+1)}\beta_{2}^{(M_{2}+2)}\cdots\beta_{2}^{(2M_{2})}\right)^{L_{2}}$,
we obtain, 
\begin{equation}
g_{M_{2}+1,\dots,2M_{2}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)+\left(\frac{\beta_{2}^{(M_{2})}}{\beta_{2}^{(M_{2}+1)}}\right)^{L_{2}}g_{M_{2},M_{2}+2,\dots,2M_{2}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)+\dots=0,\label{seq:simplified-eq}
\end{equation}
where the omitted terms are less than the order of $\left(\beta_{2}^{(M_{2})}/\beta_{2}^{(M_{2}+1)}\right)^{L_{2}}$.
Note that $\beta_{2}^{(j)},j=1,2,\dots,2M_{2}$ are dependent on $E$
and $\beta_{1}$ by $f\left(E,\beta_{1},\beta_{2}^{(j)}\right)=0$,
when $E$ is fixed, the functions $g_{j_{1}j_{2}\dots j_{M_{2}}}$
are univariate functions of $\beta_{1}$. Therefore, if $\left|\beta_{2}^{(M_{2})}\right|\neq\left|\beta_{2}^{(M_{2}+1)}\right|$,
all the terms in Eq. (\ref{seq:simplified-eq}) vanish except for
$g_{M_{2}+1,\dots,2M_{2}}\left(E,\beta_{1},\beta_{2}^{(1)},\dots,\beta_{2}^{(2M_{2})}\right)$,
so that only a finite number of solutions (that is, independent of
$L_{2}$) of $\beta_{1}$ can be solved from Eq. (\ref{seq:simplified-eq}).
Otherwise, if $\left|\beta_{2}^{(M_{2})}\right|=\left|\beta_{2}^{(M_{2}+1)}\right|$,
the first two terms in Eq. (\ref{seq:simplified-eq}) are preserved,
and the number of solutions increases with the order of $L_{2}$ when
$L_{2}\rightarrow\infty$. This is the reason why the non-Bloch waves
satisfying GBZ constraints account for a major amount of OBC eigenstates.

Given a reference energy $E$, we need to know how these solutions
are distributed on the 1D PMGBZ. Consider the relative argument $\exp\left(i\phi\right)\equiv\beta_{2}^{(M_{2}+1)}/\beta_{2}^{(M_{2})}$.
When $\phi$ is determined, the possible values of $\beta_{1}$ are
also determined. Consequently, all functions $g_{j_{1}j_{2}\dots j_{M_{2}}}$
can be viewed as local functions of $\exp\left(i\phi\right)$. In
the thermodynamic limit, Eq. (\ref{seq:simplified-eq}) reads, 
\begin{equation}
\phi=\frac{i}{L_{2}}\ln\frac{g_{M_{2},M_{2}+2,\dots,2M_{2}}\left(e^{i\phi}\right)}{g_{M_{2}+1,\dots,2M_{2}}\left(e^{i\phi}\right)}+\frac{2n\pi}{L_{2}},n=0,\pm1,\pm2,\dots.\label{seq:phi-eq}
\end{equation}
Next, we define $\tilde{g}\left(e^{i\phi}\right)\equiv i\ln\left[g_{M_{2},M_{2}+2,\dots,2M_{2}}\left(e^{i\phi}\right)/g_{M_{2}+1,\dots,2M_{2}}\left(e^{i\phi}\right)\right]$,
and consider two adjacent solutions, 
\begin{align}
\phi_{1} & =\frac{1}{L_{2}}\tilde{g}\left(e^{i\phi_{1}}\right)+\frac{2n\pi}{L_{2}},\\
\phi_{2} & =\frac{1}{L_{2}}\tilde{g}\left(e^{i\phi_{2}}\right)+\frac{2\left(n+1\right)\pi}{L_{2}},
\end{align}
then the difference of the two solutions reads, 
\begin{align}
\phi_{2}-\phi_{1} & =\frac{2\pi}{L_{2}}+\frac{1}{L_{2}}\frac{d\tilde{g}\left(\phi_{1}\right)}{d\phi}\left(\phi_{2}-\phi_{1}\right)+O\left(\frac{1}{L_{2}^{3}}\right)=\frac{2\pi}{L_{2}}+O\left(\frac{1}{L_{2}^{2}}\right),
\end{align}
where the second equation holds because the order of $\phi_{2}-\phi_{1}$
is $O(1/L_{2})$. Therefore, if the relative phase changes by $\Delta\phi$,
the number of the $\beta_{1}$-solutions reads, 
\begin{align}
N_{\text{q1D}} & =\frac{\left|\Delta\phi\right|}{\frac{2\pi}{L_{2}}+O\left(\frac{1}{L_{2}^{2}}\right)}=\frac{L_{2}}{2\pi}\left|\Delta\phi\right|+O\left(1\right),
\end{align}
which proves the relation mentioned in the main text.

\section{Topological Properties of Loop Winding Numbers and Equivalent Formulation
of $W_{\text{strip}}(E,r)$}

\label{app:winding-general-loop}


In the main text, we have introduced the strip winding number $W_{\text{strip}}\left(E,r\right)$
defined as the follows, 
\begin{equation}
W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\int_{-\pi}^{\pi}w\left(\theta_{2}\right)d\theta_{2},\label{seq:def-W-strip}
\end{equation}
where the loop winding number $w\left(\theta_{2}\right)$ is defined
as, 
\begin{equation}
w\left(\theta_{2}\right)=\frac{1}{2\pi i}\int_{-\pi}^{\pi}d\theta_{1}\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\theta_{2}}\right)\right]}{\partial\theta_{1}}.\label{seq:def-W-theta2}
\end{equation}
From a geometric perspective, $w\left(\theta_{2}\right)$ is the winding
number of the 2D ChP $f$ when $\theta_{1}$ runs over $2\pi$ and
$\theta_{2}$ keeps constant. However, in some theoretical derivations,
the exact form of $w\left(\theta_{2}\right)$ is difficult to calculate.
Instead, the winding numbers around some irregular closed loops are
rather simpler. Therefore, a more flexible version of Eq. (\ref{seq:def-W-strip})
is helpful in these cases. In this section, we will show that the
winding number $w\left(\theta_{2}\right)$ in Eq. (\ref{seq:def-W-strip})
can be substituted by $\tilde{w}\left(s\right)$ defined as, 
\begin{equation}
\tilde{w}\left(s\right)=\frac{1}{2\pi i}\int_{-\pi}^{\pi}d\theta_{1}\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\left[s+\delta_{2}\left(\theta_{1}\right)\right]}\right)\right]}{\partial\theta_{1}},\label{seq:def-tilde-w}
\end{equation}
where $\delta_{2}\left(\theta_{1}\right),\theta_{1}\in\left[-\pi,\pi\right]$
is an arbitrary function of $\theta_{1}$ satisfying $\exp\left[i\delta_{2}\left(-\pi\right)\right]=\exp\left[i\delta_{2}\left(\pi\right)\right]$.


\begin{figure}
\centering

\includegraphics{Figures/winding-general-loop-20250527}

\caption{Topological properties and equivalence of winding numbers. (a) Sketch
of the loop winding numbers as a homomorphism between fundamental
groups. (b) Relation between the loop winding numbers $w$ and $\tilde{w}$.
The red and green regions denote the regions where $\tilde{w}(\theta_{2})=w(\theta_{2})+1$
and $\tilde{w}(\theta_{2})=w(\theta_{2})-1$, respectively. (c) Simplification
of the winding number around the loop that winds around $\theta_{2}$-axis.
(d) Homology transformation of an $n$-fold winding loop (blue curve
in i) to $n$ congruent one-fold winding loops (blue solid curve and
purple dashed curve in ii).}
\label{sfig:winding-general-loop} 
\end{figure}

Before discussing the winding number $\tilde{w}\left(s\right)$, we
first give an overview of the topological properties of the winding
number on the toric base space $X(E,r)$. As narrated in the main
text, the toric base space is defined as Eq. (\ref{seq:toric-base-space}),
\begin{equation}
X(E,r)=\left\{ \left(re^{i\theta_{1}},\rho\left(\theta_{1}\right)e^{i\theta_{2}}\right)\in\mathbb{C}^{2}|\theta_{1},\theta_{2}\in\left[-\pi,\pi\right]\right\} ,\label{seq:toric-base-space}
\end{equation}
which is homeomorphic to the 2D torus $\mathbb{T}^{2}$. For a closed
loop $\gamma$ on the space $X$ defined by the parametric equations
$\theta_{1}=\theta_{1}\left(t\right),\theta_{2}=\theta_{2}\left(t\right),t\in\left[0,1\right]$
satisfying $\exp\left[i\theta_{1}\left(0\right)\right]=\exp\left[i\theta_{1}\left(1\right)\right]$
and $\exp\left[i\theta_{2}\left(0\right)\right]=\exp\left[i\theta_{2}\left(1\right)\right]$,
we can define the general loop winding number $w_{\text{loop}}$ defined
as, 
\begin{equation}
w_{\text{loop}}\left(\gamma\right)\equiv\frac{1}{2\pi i}\int_{0}^{1}dt\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}(t)},\rho\left(\theta_{1}(t)\right)e^{i\theta_{2}(t)}\right)\right]}{\partial t}.\label{seq:w-loop}
\end{equation}
The geometric picture of $w_{\text{loop}}$ is illustrated in Fig.
\ref{sfig:winding-general-loop}(a). Because the map $f\left(E,\cdot,\cdot\right):X(E,r)\rightarrow\mathbb{C}$
is a continuous map, closed loops {[}loops $\gamma_{1}$ and $\gamma_{2}$
in Fig. \ref{sfig:winding-general-loop}(a) {]} on $X(E,r)$ are mapped
to closed loops on $\mathbb{C}$ {[}purple and orange closed loops
in the complex plane of Fig. \ref{sfig:winding-general-loop}(a) {]}.
When the image of $\gamma$ under $f\left(E,\cdot,\cdot\right)$ does
not pass $0$, $w_{\text{loop}}\left(\gamma\right)$ is well defined,
and equal to the number of laps that the image of $\gamma$ winds
around $0$. Because the preimage of $0$ under $f\left(E,\cdot,\cdot\right)$
is $\text{PMGBZ}\left(E\right)\cap X(E,r)$, we consider the restriction
of $f\left(E,\cdot,\cdot\right)$ on $X(E,r)\backslash\text{PMGBZ}\left(E\right)$,
that is, $f\left(E,\cdot,\cdot\right):X(E,r)\backslash\text{PMGBZ}\left(E\right)\rightarrow\mathbb{C}\backslash\left\{ 0\right\} $.
When $\text{PMGBZ}\left(E\right)\cap X(E,r)$ is a finite set, $X(E,r)\backslash\text{PMGBZ}\left(E\right)$
is pathwise connected, so that $f\left(E,\cdot,\cdot\right)$ induces
a homomorphism $f_{*}\left(E,\cdot,\cdot\right):\pi_{1}\left(X(E,r)\backslash\text{PMGBZ}\left(E\right)\right)\rightarrow\pi_{1}\left(\mathbb{C}\backslash\left\{ 0\right\} \right)=\pi_{1}\left(\mathbb{T}^{1}\right)=\mathbb{Z}$,
which maps the homotopy class $\left[\gamma\right]$ to $w_{\text{loop}}\left(\gamma\right)$
\cite{hatcherAlgebraicTopology2001}. In other words, the map $w_{\text{loop}}:\left\{ \text{Closed loops on }X(E,r)\backslash\text{PMGBZ}\left(E\right)\right\} \rightarrow\mathbb{Z}$
satisfies $w_{\text{loop}}\left(\gamma_{1}\circ\gamma_{2}\right)=w_{\text{loop}}\left(\gamma_{1}\right)+w_{\text{loop}}\left(\gamma_{2}\right)$,
where $\circ$ is the product of the fundamental group $\pi_{1}\left(X(E,r)\backslash\text{PMGBZ}\left(E\right)\right)$.
Furthermore, because $\mathbb{Z}$ is abelian, $w_{\text{loop}}$
vanishes on the commutators, that is, $w_{\text{loop}}\left(\gamma_{1}\circ\gamma_{2}\circ\gamma_{1}^{-1}\circ\gamma_{2}^{-1}\right)=0,\forall\gamma_{1},\gamma_{2}\in X(E,r)\backslash\text{PMGBZ}(E)$,
so that $w_{\text{loop}}$ is also a homomorphism from the homology
group $H_{1}\left(X(E,r)\backslash\text{mGBZ}\left(E\right)\right)$
\cite{hatcherAlgebraicTopology2001} to the integer group $\mathbb{Z}$.

The topological nature of $w_{\text{loop}}$ allows us to perform
algebraic operations on the winding numbers of the closed loops. As
an example, the difference $\gamma_{2}-\gamma_{1}$ (in the sense
of homology) in Fig. \ref{sfig:winding-general-loop}(a) is homologous
to an infinitesimal loop around the PMGBZ point, shown as the red
arrowed loop in Fig. \ref{sfig:winding-general-loop}(a). As a result,
the difference $w_{\text{loop}}\left(\gamma_{2}\right)-w_{\text{loop}}\left(\gamma_{1}\right)$
is equal to the winding number around the infinitesimal loop, which
is defined as the topological charge of the PMGBZ point.

Based on the discussion above, we will show that $W_{\text{strip}}$
can also be calculated by the following integral, 
\begin{equation}
W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\int_{-\pi}^{\pi}\tilde{w}\left(s\right)ds.\label{seq:W-strip-new}
\end{equation}
As defined in Eq. (\ref{seq:def-tilde-w}), for each constant value
$s$, $\tilde{w}\left(s\right)$ equals the winding number of the
loop defined by $\theta_{2}\left(\theta_{1}\right)=\delta_{2}\left(\theta_{1}\right)+s$
and $\theta_{1}=2\pi t,t\in\left[0,1\right]$. We first consider the
case where $\delta_{2}\left(-\pi\right)=\delta_{2}\left(\pi\right)$,
that is, the loop does not wind around the $\theta_{2}$-axis. Without
loss of generality, we suppose $\delta_{2}(-\pi)=\delta_{2}(\pi)=0$.
As shown in Fig. \ref{sfig:winding-general-loop}(b), for each constant
value $s$, the difference between $\tilde{w}\left(s\right)$ and
$w\left(s\right)$ equals the total topological charges enclosed by
the loop of $\tilde{w}\left(s\right)$ {[}orange solid line in Fig.
\ref{sfig:winding-general-loop}(b){]} and the reverse of the loop
of $w\left(s\right)$ {[}orange dashed line in Fig. \ref{sfig:winding-general-loop}(b){]}.
For a pair of PMGBZ points, as illustrated by the light red and light
green regions in Fig. \ref{sfig:winding-general-loop}(b), the contributions
of the positive charge and the negative charge to $\int_{-\pi}^{\pi}\left[\tilde{w}\left(s\right)-w\left(s\right)\right]ds$
cancel with each other, so that the integral of $\tilde{w}\left(s\right)$
equals the integral of $w\left(s\right)$, which proves Eq. (\ref{seq:W-strip-new}).

Next, we consider the general case where the loop defined by $\theta_{2}=\delta_{2}\left(\theta_{1}\right)$
is allowed to wind around the $\theta_{2}$-axis, shown as the blue
solid line in Fig. \ref{sfig:winding-general-loop}(c). As illustrated
in the figure, the blue loop is homotopic to the purple dashed loop,
and the purple dashed loop is homologous to the sum of the orange
dotted loop, which satisfies the condition $\delta_{2}\left(-\pi\right)=\delta_{2}\left(\pi\right)$,
and the gray loop, around which the winding number of ChP vanishes.
Therefore, for each winding loop (solid blue line), the winding number
around the loop equals the winding number of a special loop satisfying
$\delta_{2}\left(-\pi\right)=\delta_{2}\left(\pi\right)$ (dotted
orange line), which is the case shown in Fig. \ref{sfig:winding-general-loop}(b).

Furthermore, in some cases, the form of the ChP is complicated, but
the product, 
\begin{equation}
g\left(\theta_{1},\theta_{2}\right)=\prod_{m=1}^{n}f\left(E,re^{i\theta_{1}},\rho(\theta_{1})e^{i\left(\theta_{2}+2m\pi/n\right)}\right),
\end{equation}
has a simple form, such as the $[11]$-SGBZ in 2D HN model discussed
in the main text. Then, $W_{\text{strip}}$ can be calculated by,
\begin{equation}
W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\int_{-\pi/n}^{\pi/n}\tilde{w}_{g}\left(s\right)ds,\label{eq:W-strip-g}
\end{equation}
where, 
\begin{equation}
\tilde{w}_{g}\left(s\right)=\int_{-\pi}^{\pi}\frac{d\theta_{1}}{2\pi i}\frac{\partial\ln\left[g\left(\theta_{1},s+\delta_{2}(\theta_{1})/n\right)\right]}{\partial\theta_{1}},\label{seq:w-g}
\end{equation}
is the winding number of $g$ on the loop $s+\delta_{2}(\theta_{1})/n$.
Here, we only require $\delta_{2}(\theta_{1})$ rather than $\delta_{2}(\theta_{1})/n$
to satisfy the periodicity condition $\exp[i\delta_{2}(-\pi)]=\exp[i\delta_{2}(\pi)]$.
Now, we will prove Eq. (\ref{seq:w-g}).

First, $g(\theta_{1},s+\delta_{2}(\theta_{1})/n)$ is a periodic function
of $\theta_{1}$. When $\theta_{1}$ increases by $2\pi$, $\delta_{2}(\theta_{1})/n$
increases by an integer multiple of $2\pi/n$. Supposing $\delta_{2}(\pi)/n=\delta_{2}(-\pi)/n+2\pi m_{0}/n\mod{2\pi}$,
we get, 
\begin{align}
g\left(\pi,s+\delta_{2}(\pi)/n\right) & =\prod_{m=1}^{n}f\left(E,re^{i\pi},\rho(\pi)e^{i\left(s+\delta_{2}(\pi)/n+2m\pi/n\right)}\right),\nonumber \\
 & =\prod_{m=1}^{n}f\left(E,re^{-i\pi},\rho(-\pi)e^{i\left(s+\delta_{2}(-\pi)/n+2(m+m_{0})\pi/n\right)}\right),\nonumber \\
 & =g\left(-\pi,s+\delta_{2}(-\pi)/n\right).
\end{align}
Therefore, the number $\tilde{w}_{g}(s)$ defined by Eq. (\ref{seq:w-g})
is a valid winding number.

Next, we consider the relation between $\tilde{w}_{g}(s)$ and the
winding number of the ChP. When $\delta_{2}(-\pi)/n=\delta_{2}(\pi)/n\mod{2\pi}$,
the curve $\theta_{2}=s+\delta_{2}(\theta_{1})/n$ itself forms a
closed loop on $X(E,r)$. Therefore, the following relation holds,
\begin{align}
\tilde{w}_{g}\left(s\right) & =\sum_{m=1}^{n}\int_{-\pi}^{\pi}\frac{d\theta_{1}}{2\pi i}\frac{\partial\ln\left[f\left(E,re^{i\theta_{1}},\rho(\theta_{1})e^{i\left(s+\delta_{2}(\theta_{1})/n+2m\pi/n\right)}\right)\right]}{\partial\theta_{1}},\nonumber \\
 & =\sum_{m=1}^{n}\tilde{w}\left(s+2m\pi/n\right),\label{seq:relation-winding-numbers}
\end{align}
where $\tilde{w}(s)$ here is defined by the loop $\delta_{2}^{\prime}(\theta_{1})\equiv\delta_{2}(\theta_{1})/n$.
By Eq. (\ref{seq:W-strip-new}), the integral of Eq. (\ref{eq:W-strip-g})
reads, 
\begin{align}
\frac{1}{2\pi}\int_{-\pi/n}^{\pi/n}\tilde{w}_{g}\left(s\right)ds= & \sum_{m=1}^{n}\frac{1}{2\pi}\int_{-\pi/n}^{\pi/n}\tilde{w}\left(s+2m\pi/n\right)ds,\nonumber \\
= & \frac{1}{2\pi}\int_{\pi/n}^{2\pi+\pi/n}\tilde{w}\left(s\right)ds,
\end{align}
which equals $W_{\text{strip}}(E,r)$.

When $\delta_{2}(-\pi)/n=\delta_{2}(\pi)/n+2\pi m_{0}/n\mod{2\pi},m_{0}\neq0$,
the curve $\theta_{2}=s+\delta_{2}(\theta_{1})/n,\theta_{1}\in[-\pi,\pi]$
is not closed, so that the winding number $\tilde{w}(s)$ in Eq. (\ref{seq:relation-winding-numbers})
is ill-defined. However, by homology transformations, $\tilde{w}_{g}(s)$
can be transformed into the sum of the winding numbers around $n$
congruent winding loops with a phase shift of $2\pi/n$ along the
$\theta_{2}$-axis. Taking $n=2,m_{0}=1$ as an example, as shown
in panel i of Fig. \ref{sfig:winding-general-loop}(d), the curve
$\delta_{2}(\theta_{1})/2$ is braided with $\delta_{2}(\theta_{1})/2+\pi$
(black dotted curves), forming a twofold winding loop around the $\theta_{1}$-axis
(solid blue curves). The winding number $\tilde{w}_{g}(s)$ is equal
to the winding number of the ChP around the twofold loop. By virtue
of the homology invariance of the winding number, a winding loop parallel
to the $\theta_{2}$-axis can be added to the two-fold loop without
changing the winding number, shown as the gray horizontal line in
panel i of Fig. \ref{sfig:winding-general-loop} (d). Then, the sum
of the two-fold loop and the horizontal loop is homologous to the
sum of two one-fold loops, shown as the blue solid curve and the purple
dashed curve in panel ii of Fig. \ref{sfig:winding-general-loop}(d).
Thus, we return to the case with $m_{0}=0$. In general, with the
same method, an arbitrary $n$-fold loop can be split into $n$ one-fold
loops by homology, so Eq. (\ref{eq:W-strip-g}) holds for arbitrary
$n$ and $m_{0}$.

\section{Sign of Topological Charge and Change of Relative Argument}

\label{app:topo-charge-change-phi}

In the main text, we show that the increment of strip winding number
is related to the change of relative arguments by, 
\begin{equation}
W_{\text{strip}}\left(E,r^{\prime}\right)-W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\sum_{j=1}^{k}\left(-1\right)^{\tau_{j}}\left(\phi_{j}^{\prime}-\phi_{j}\right),\label{seq:W-strip}
\end{equation}
where $\phi_{j}^{\prime}$ and $\phi_{j}$ are the relative arguments
of the $j$-th PMGBZ pair at radius $r$ and $r^{\prime}$, and $\tau_{j}=0,1$
is related to the topological charge of the PMGBZ pair. If the base
point of the relative phase $\phi_{j}$ takes a positive charge, $\tau_{j}=0$,
and vice versa. In this section, we will show that the sign of $\phi_{j}^{\prime}-\phi_{j}$
is related to the topological charge when $r^{\prime}>r$.

For an PMGBZ pair consisting of two PMGBZ points $p_{\mathrm{a}}=\left(\mathcal{B}_{1},\mathcal{B}_{2}\right)$
and $p_{\mathrm{b}}=\left(\mathcal{B}_{1},\mathcal{B}_{2}e^{i\phi}\right)$,
as schematically illustrated in Fig. \ref{sfig:topo-charge-rel-arg}(a),
there must exist two curves of $\beta_{2}$-zeros of 2D ChP in the
two neighborhoods of $p_{\mathrm{a}}$ and $p_{\mathrm{b}}$, shown
as the curves $\beta_{2}^{(\mathrm{a})}$ and $\beta_{2}^{(\mathrm{b})}$
in Fig. \ref{sfig:topo-charge-rel-arg}(a). It is noted that the ordering
$|\beta_{2}^{(1)}|\leq|\beta_{2}^{(2)}|\leq\cdots\leq|\beta_{2}^{(M_{2}+N_{2})}|$
used in the main text will lead to the discontinuity of $\beta_{2}^{(j)}$
as a function of $\theta_{1}$. For example, the $M_{2}$-th solutions
{[}blue and brown dashed curves in Fig. \ref{sfig:topo-charge-rel-arg}(a){]}
jump into $\left(M_{2}+1\right)$-th solutions {[}blue and brown solid
curves in Fig. \ref{sfig:topo-charge-rel-arg}(a){]} when passing
an PMGBZ point. Therefore, in this section, we use the superscript
$(\mathrm{a})$ or $(\mathrm{b})$ to label the continuous $\beta_{2}$-zero
curve that passes the PMGBZ point $p_{\mathrm{a}}$ or $p_{\mathrm{b}}$,
instead of the indices ordered by absolute values.

\begin{figure}
\centering

\includegraphics{Figures/supp-derivative-zeros-20250120}

\caption{Sketch for the relation between topological charge and relative phase.
(a) 3D view of a PMGBZ pair and the $\beta_{2}$-zeros in the $\theta_{1}$-$\beta_{2}$
space. The gray surface is the toric base manifold $X(E,r)$, $p_{\mathrm{a}}$
and $p_{\mathrm{b}}$ form a PMGBZ pair, and $\beta_{2}^{(\mathrm{a})},\beta_{2}^{(\mathrm{b})}$
are the $\beta_{2}$-zeros of $f(E,re^{i\theta_{1}},\beta_{2}^{(j)})$
in the neighborhoods of $p_{\mathrm{a}}$ and $p_{\mathrm{b}}$, respectively.
(b) Absolute value of $\beta_{2}^{(\mathrm{a})}$ and $\beta_{2}^{(\mathrm{b})}$
in the neighborhood of a PMGBZ pair.}
\label{sfig:topo-charge-rel-arg} 
\end{figure}


The topological charge of a PMGBZ point is determined by the direction
in which the $\beta_{2}$-zero curve passes through the base space.
As illustrated in Fig. \ref{sfig:topo-charge-rel-arg}(b), because
the radius function $\rho\left(\theta_{1}\right)$ must be sandwiched
between $|\beta_{2}^{(\mathrm{a})}|$ and $|\beta_{2}^{(\mathrm{b})}|$,
the direction of the $\beta_{2}$-zero curves is determined by the
derivative of $\left|\beta_{2}^{(\mathrm{a})}\right|$ and $\left|\beta_{2}^{(\mathrm{b})}\right|$
against $\theta_{1}$ at the PMGBZ points. If $\frac{d\left|\beta_{2}^{(\mathrm{a})}\right|^{2}}{d\theta_{1}}|_{p_{\mathrm{a}}}>\frac{d\left|\beta_{2}^{(\mathrm{b})}\right|^{2}}{d\theta_{1}}|_{p_{\mathrm{b}}}$,
$p_{\mathrm{a}}$ takes the positive charge and $p_{\mathrm{b}}$
the negative charge, and vice versa. In general, the derivative of
a $\beta_{2}$-zero is derived from the implicit equation $f\left(E,re^{i\theta_{1}},\beta_{2}^{(u)}\right)=0$,
i.e., 
\begin{equation}
\frac{d\beta_{2}^{(u)}}{d\theta_{1}}=-i\frac{\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{u}}re^{i\theta_{1}}}{\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{u}}},\label{seq:derivative-beta-2}
\end{equation}
where $u=\mathrm{a},\mathrm{b}$. Then, the derivative of the absolute
value reads, 
\begin{align}
\frac{d\left|\beta_{2}^{(u)}\right|^{2}}{d\theta_{1}} & =2\text{Re}\left(\beta_{2}^{(u)*}\frac{d\beta_{2}^{(u)}}{d\theta_{1}}\right)=2\text{Im}\left(\frac{\beta_{2}^{(u)*}re^{i\theta_{1}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{u}}}{\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{u}}}\right).
\end{align}
For simplicity, we denote 
\begin{equation}
K\left(E,\beta_{1},\beta_{2}\right)=\frac{\beta_{2}^{*}\beta_{1}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{\left(E,\beta_{1},\beta_{2}\right)}}{\left.\frac{\partial f}{\partial\beta_{2}}\right|_{\left(E,\beta_{1},\beta_{2}\right)}},\label{eq:K-fun}
\end{equation}
then we have $\frac{d\left|\beta_{2}^{(\mathrm{a})}\right|^{2}}{d\theta_{1}}|_{p_{\mathrm{a}}}=2\text{Im}\left[K\left(E,\mathcal{B}_{1},\mathcal{B}_{2}\right)\right]$
and $\frac{d\left|\beta_{2}^{(\mathrm{b})}\right|^{2}}{d\theta_{1}}|_{p_{\mathrm{b}}}=2\text{Im}\left[K\left(E,\mathcal{B}_{1},\mathcal{B}_{2}e^{i\phi}\right)\right]$.
By comparing the derivatives, the topological charges of $p_{\mathrm{a}}$
and $p_{\mathrm{b}}$ are determined.

Next, we will show the relation between the function $K$ defined
above and the change in the relative phase $\phi$ when the modulus
of $\mathcal{B}_{1}$ increases, i.e., the sign of $\partial\left|\mathcal{B}_{1}\right|^{2}/\partial\phi$
when $E$ is kept constant. The derivative can be calculated from
the auxiliary GBZ equations, 
\begin{align}
f\left(E,\mathcal{B}_{1},\mathcal{B}_{2}\right) & =0,\label{seq:aGBZ-1}\\
f\left(E,\mathcal{B}_{1},\mathcal{B}_{2}e^{i\phi}\right) & =0,\label{seq:aGBZ-2}
\end{align}
as implicit functions of $\mathcal{B}_{1}\left(E,\phi\right)$ and
$\mathcal{B}_{2}\left(E,\phi\right)$. By direct calculations, the
derivative reads, 
\begin{equation}
\frac{\partial\mathcal{B}_{1}}{\partial\phi}=\frac{-ie^{i\phi}\mathcal{B}_{2}\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{a}}}\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{b}}}}{\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{a}}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{\mathrm{b}}}-e^{i\phi}\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{b}}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{\mathrm{a}}}},\label{eq:diff-mGBZ-beta1}
\end{equation}
and the derivative of $\left|\mathcal{B}_{1}\right|^{2}$ reads, 
\begin{align}
\frac{\partial\left|\mathcal{B}_{1}\right|^{2}}{\partial\phi} & =2C\text{Im}\left[K\left(E,\mathcal{B}_{1},\mathcal{B}_{2}\right)-K\left(E,\mathcal{B}_{1},\mathcal{B}_{2}e^{i\phi}\right)\right],\nonumber \\
 & =C\left(\left.\frac{d|\beta_{2}^{(\mathrm{a})}|^{2}}{d\theta_{1}}\right|_{p_{\mathrm{a}}}-\left.\frac{d|\beta_{2}^{(\mathrm{b})}|^{2}}{d\theta_{1}}\right|_{p_{\mathrm{b}}}\right),\label{seq:relation-dnorm-K-fun}
\end{align}
where, 
\begin{equation}
C=\frac{\left|\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{b}}}\right|^{2}\cdot\left|\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{a}}}\right|^{2}}{\left|\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{a}}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{\mathrm{b}}}-e^{i\phi}\left.\frac{\partial f}{\partial\beta_{2}}\right|_{p_{\mathrm{b}}}\left.\frac{\partial f}{\partial\beta_{1}}\right|_{p_{\mathrm{a}}}\right|^{2}},
\end{equation}
is a positive value. Therefore, the sign of $\partial\left|\mathcal{B}_{1}\right|^{2}/\partial\phi$
is the same as the topological charge of $p_{\mathrm{a}}$, and opposite
to the topological charge of $p_{\mathrm{b}}$.

For the strip winding number $W_{\text{strip}}\left(E,r\right)$,
the sign $\left(-1\right)^{\tau_{j}}$ equals $+1$ when $\phi_{j}$
starts from a positive charge, and equals $-1$ when starting from
a negative charge. According to the discussion above, $\left(-1\right)^{\tau_{j}}$
equals the sign of $\partial\left|\mathcal{B}_{1}\right|^{2}/\partial\phi$.
That is, when the radius $r$ increases to some $r^{\prime}>r$, the
sign of $\phi_{j}^{\prime}-\phi_{j}$ equals $\left(-1\right)^{\tau_{j}}$.
As a result, Eq. (\ref{seq:W-strip}) is simplified into, 
\begin{equation}
W_{\text{strip}}\left(E,r^{\prime}\right)-W_{\text{strip}}\left(E,r\right)=\frac{1}{2\pi}\sum_{j=1}^{k}\left|\phi_{j}^{\prime}-\phi_{j}\right|,
\end{equation}
which proves the claim in the main text.

\section{Relation of the Degrees in Quasi-1D and 2D Characteristic Polynomials}

\label{app:relation-degrees}

In the main text, we state without proof that the lowest and highest
degrees of $\beta_{1}$ in the quasi-1D ChP $F(E,\beta_{1})$ and
2D ChP $f(E,\beta_{1},\beta_{2})$ are related to each other by, 
\begin{align}
M & =M_{1}L_{2}+O(1),\label{seq:relation-M}\\
N & =N_{1}L_{2}+O(1).\label{seq:relation-N}
\end{align}
In this section, we will prove the two relations for a general 2D
lattice.

First, we expand the 2D non-Bloch Hamiltonian $h(\beta_{1},\beta_{2})$
as the polynomial of $\beta_{2}$, i.e., 
\begin{equation}
h\left(\beta_{1},\beta_{2}\right)=\sum_{j=-M_{2}}^{N_{2}}h^{(j)}\left(\beta_{1}\right)\beta_{2}^{j},\label{seq:h-in-polynomial}
\end{equation}
where $h^{(j)}(\beta_{1})$ is an $m\times m$ matrix. Then, the hybrid
non-Bloch Hamiltonian can be expanded into the block Toeplitz form
under the single particle basis, i.e., 
\begin{equation}
\mathcal{H}\left(\beta_{1}\right)=\underbrace{\left(\begin{matrix}h^{(0)} & h^{(-1)} & \cdots & h^{(-M_{2})}\\
h^{(1)} & h^{(0)} & h^{(-1)} & \cdots & \ddots\\
\vdots & h^{(1)} & h^{(0)} & h^{(-1)} & \cdots & h^{(-M_{2})}\\
h^{(N_{2})} & \vdots & h^{(1)} & h^{(0)} & \ddots & \vdots\\
 & \ddots & \vdots & \ddots & \ddots & h^{(-1)}\\
 &  & h^{(N_{2})} & \cdots & h^{(1)} & h^{(0)}
\end{matrix}\right)}_{L_{2}\text{ blocks}}.\label{seq:quasi-1D-Toeplitz}
\end{equation}

We first consider the lowest degrees. The case of the highest degrees
can be proved by the same reasoning. It is noted that the matrix elements
of $h^{(j)}(\beta_{1})$ are all Laurent polynomials of $\beta_{1}$.
Assume that the lowest degree of $\beta_{1}$ in $h_{\mu\nu}(\beta_{1},\beta_{2})$
is $-M_{\mu\nu}$, and the degree of $\beta_{2}$ in the term with
lowest degree of $\beta_{1}$ is $t_{\mu\nu}$, then, the elements
of the ChP $f(E,\beta_{1},\beta_{2})$ can be expressed as the leading
term $a_{ij}\beta_{1}^{-M_{ij}}\beta_{2}^{t_{ij}}$ plus a remainder
in which the degrees of $\beta_{1}$ are greater than $-M_{ij}$,
like Eq. (\ref{seq:expansion-f}), 
\begin{align}
f\left(E,\beta_{1},\beta_{2}\right) & =\det\left[E-h\left(\beta_{1},\beta_{2}\right)\right],\nonumber \\
 & =\begin{vmatrix}a_{11}\beta_{1}^{-M_{11}}\beta_{2}^{t_{11}}+\cdots & a_{12}\beta_{1}^{-M_{12}}\beta_{2}^{t_{12}}+\cdots & \cdots & a_{1m}\beta_{1}^{-M_{1m}}\beta_{2}^{t_{1m}}+\cdots\\
a_{21}\beta_{1}^{-M_{21}}\beta_{2}^{t_{21}}+\cdots & a_{22}\beta_{1}^{-M_{22}}\beta_{2}^{t_{22}}+\cdots &  & \vdots\\
\vdots &  & \ddots\\
a_{m1}\beta_{1}^{-M_{m1}}\beta_{2}^{t_{m1}}+\cdots & \cdots &  & a_{mm}\beta_{1}^{-M_{mm}}\beta_{2}^{t_{mm}}+\cdots
\end{vmatrix},\label{seq:expansion-f}
\end{align}
where $a_{ij}\in\mathbb{C}$ are the coefficients. In the determinant
of Eq. (\ref{seq:expansion-f}), only the leading terms contribute
to the lowest degree of $f(E,\beta_{1},\beta_{2})$. Therefore, the
term in $f(E,\beta_{1},\beta_{2})$ with lowest degree in $\beta_{1}$
can be construct by the following steps: 
\begin{enumerate}
\item Find $i_{\max},j_{\max}$ that maximizes $M_{ij}$, i.e. $M_{i_{\max}j_{\max}}=\max_{i,j}\left\{ M_{ij}\right\} $. 
\item Get the $(i_{\max},j_{\max})$-cofactor of $f\left(E,\beta_{1},\beta_{2}\right)$,
denoted as $f_{1}(E,\beta_{1},\beta_{2})$, then factorize $f(E,\beta_{1},\beta_{2})$
as, 
\begin{equation}
f\left(E,\beta_{1},\beta_{2}\right)=a_{i_{\max}j_{\max}}\beta_{1}^{-M_{i_{\max}j_{\max}}}\beta_{2}^{t_{i_{\max}j_{\max}}}f_{1}\left(E,\beta_{1},\beta_{2}\right)+\dots
\end{equation}
\item Substitute the determinant $f$ by the cofactor $f_{1}$, and repeat
the two steps above to get $f_{2}$, $f_{3}$, $\dots$, until the
order of the cofactor is reduced to $1$. 
\end{enumerate}
After this procedure, we get a series of indices $(i_{\max}^{(1)},j_{\max}^{(1)}),(i_{\max}^{(2)},j_{\max}^{(2)}),\dots,(i_{\max}^{(m)},j_{\max}^{(m)})$.
According to our construction, $i_{\max}^{(1)},i_{\max}^{(2)},\dots,i_{\max}^{(m)}$
and $j_{\max}^{(1)},j_{\max}^{(2)},\dots,j_{\max}^{(m)}$ are permutations
of $1,2,\dots,m$. Therefore, we can sort $M_{i_{\max}^{(1)}j_{\max}^{(1)}},M_{i_{\max}^{(2)}j_{\max}^{(2)}},\dots,M_{i_{\max}^{(m)}j_{\max}^{(m)}}$
to $M_{1\nu_{1}},M_{2\nu_{2}},\dots,M_{m\nu_{m}}$, and ensure that
$\nu_{1},\nu_{2},\dots,\nu_{m}$ is a permutation of $1,2,\dots,m$.
Then, the minimum negative degree of $f(E,\beta_{1},\beta_{2})$ reads,
\begin{equation}
M_{1}=\sum_{j=1}^{m}M_{j\nu_{j}}.
\end{equation}

Back to the hybrid non-Bloch Hamiltonian $\mathcal{H}(\beta_{1})$
of Eq. (\ref{seq:quasi-1D-Toeplitz}), according to Eq. (\ref{seq:h-in-polynomial}),
the term $a_{j\nu_{j}}\beta_{1}^{-M_{j\nu_{j}}}$ is available in
the matrix element $h_{j\nu_{j}}^{(t_{j\nu_{j}})}(\beta_{1})$. Therefore,
in each row of the blocks in Eq. (\ref{seq:quasi-1D-Toeplitz}), we
can always pick the terms $a_{1\nu_{1}}\beta_{1}^{-M_{1\nu_{1}}},a_{2\nu_{2}}\beta_{2}^{-M_{2\nu_{2}}},\dots,a_{m\nu_{m}}\beta_{2}^{-M_{m\nu_{m}}}$,
except for the first $N_{2}$ or last $M_{2}$ rows. For different
rows of the blocks, the column indices of the selected terms do not
coincide with each other, in that $\nu_{1},\nu_{2},\dots,\nu_{m}$
is a permutation of $1,2,\dots,m$. Next, considering the determinant
$F(E,\beta_{1})\equiv\det[E-\mathcal{H}(\beta_{1})]$, except for
the first $N_{2}$ and last $M_{2}$ rows of the blocks, the lowest
degree of $\beta_{1}$ contributed from each row reads, 
\[
\beta_{1}^{-M_{1\nu_{1}}}\beta_{1}^{-M_{2\nu_{2}}}\cdots\beta_{1}^{-M_{m\nu_{m}}}=\beta_{1}^{-M_{1}},
\]
and the total contribution of the $L_{2}-M_{2}-N_{2}$ rows is, 
\[
\beta_{1}^{-M_{1}\left(L_{2}-M_{2}-N_{2}\right)}=\beta_{1}^{-M_{1}L_{2}+O(1)}.
\]
Because the terms of $\beta_{1}$ in the first $N_{2}$ and last $M_{2}$
rows can only change the total degree of $\beta_{1}$ by a finite
amount (i.e., independent of $L_{2}$), the lowest degree of $F(E,\beta_{1})$
satisfies Eq. (\ref{seq:relation-M}). By the same reasoning, we can
also prove Eq. (\ref{seq:relation-N}).

\section{Numerical Calculation of Quasi-1D Winding Number and Strip Winding
Number}

\label{app:num-W-strip}


In the main text, we have derived the relation between the quasi-1D
winding number $W(E,r)$ and the strip winding number $W_{\text{strip}}(E,r)$,
which reads $W_{\text{strip}}(E,r)=W(E,r)/L_{2}+O(1/L_{2})$. In order
to illustrate the relation, we numerically calculate $W(E,r)$ and
$W_{\text{strip}}(E,r)$ in two specific models: the 2D HN model and
the non-Hermitian Haldane model, and compare the curves of $W(E,r)/L_{2}$
and $W_{\text{strip}}(E,r)$ for randomly selected reference energies.


\begin{figure}
\centering

\includegraphics{Figures/winding-number-numerical-20250212}

\caption{Comparison of $W(E,r)$ and $W_{\text{strip}}(E,r)$ in 2D HN model
and non-Hermitian Haldane model. (a) Sketch of 2D HN model, where
the parameters are $J_{x1}=1+i$, $J_{x2}=1.5+1.2i$, $J_{y1}=-1+i$
and $J_{y2}=-1.2-0.5i$. The major (minor) axis is marked by the lattice
vector $\mathbf{a}_{1}$ ($\mathbf{a}_{2}$). (b) Numerical calculation
of $W(E,r)/L_{2}$ (blue solid lines) and $W_{\text{strip}}(E,r)$
(orange dashed lines) of 2D HN model with major axis $\mathbf{a}_{1}$.
The width is $L_{2}=40$, and the reference energies $E$ are (i)
$1.00296-0.21641i$, (ii) $1.55832+0.91741i$, (iii) $-2.57608+1.13451i$
and (iv) $-1.57168+0.02125i$. (c) Sketch of non-Hermitian Haldane
model, where the parameters are $t_{1}=0.70502$, $t_{2}=-1.32760$,
$\gamma=2.15618$, $\varphi=0.05877$ and $m=-0.64569$. The unit
cell is marked by the purple ellipse, and the major (minor) axis is
marked by the lattice vector $\mathbf{a}_{1}$ ($\mathbf{a}_{2}$).
(d) Numerical calculation of $W(E,r)/L_{2}$ (blue solid lines) and
$W_{\text{strip}}(E,r)$ (orange dashed lines) of the non-Hermitian
Haldane model with major axis $\mathbf{a}_{1}$. The width is $L_{2}=20$
($40$ sites in the supercell), and the reference energies $E$ are
(i) $0.01162+0.68736i$, (ii) $0.54100-1.99811i$, (iii) $0.59046+1.89903i$
and (iv) $0.10591-0.34594i$.}
\label{sfig:numerical-comparison-winding} 
\end{figure}

Figure \ref{sfig:numerical-comparison-winding}(a) shows the 2D HN
model with complex coupling coefficients, the same as the example
in the main text. In the numerical calculation, the coefficients are
the same as the numerical examples of the main text, which are $J_{x1}=1+i$,
$J_{x2}=1.5+1.2i$, $J_{y1}=-1+i$ and $J_{y2}=-1.2-0.5i$. The lattice
vector $\mathbf{a}_{1}=(1,1)$ is selected as the main axis and $\mathbf{a}_{2}=(0,1)$
as the minor axis. To calculate the quasi-1D winding number, the supercell
is selected by selecting successive $L_{2}$ unit cells along the
minor axis $\mathbf{a}_{2}$. Figure \ref{sfig:numerical-comparison-winding}(b)
illustrates the numerical result of $W_{\text{strip}}(E,r)$ and $W(E,r)/L_{2}$
of the 2D HN model with randomly generated reference energies. In
the numerical calculation, $L_{2}$ is set to $40$, and the reference
energies are (i) $1.00296-0.21641i$, (ii) $1.55832+0.91741i$, (iii)
$-2.57608+1.13451i$ and (iv) $-1.57168+0.02125i$. As shown in Fig.
\ref{sfig:numerical-comparison-winding}(b), the curves of $W_{\text{strip}}(E,r)$
against $r$ (orange dashed curves) are piecewise smooth curves, while
the curves of $W(E,r)/L_{2}$ (blue solid curves) show platforms because
$L_{2}$ is finite. For all four samples in Fig. \ref{sfig:numerical-comparison-winding}(b),
$W(E,r)/L_{2}$ fits well with $W_{\text{strip}}(E,r)$.

To show the generality of our conclusion, we also compare $W(E,r)$
and $W_{\text{strip}}(E,r)$ in a non-Hermitian version of the Haldane
model. As shown in Fig. \ref{sfig:numerical-comparison-winding}(c),
the nearest coupling is $t_{1}\in\mathbb{R}$, and the next-nearest
coupling has a nonreciprocal phase, i.e., the coupling coefficient
is $t_{2}e^{i\left(\varphi+\gamma\right)}$ along the direction of
the arrow, and $t_{2}e^{-i\left(\varphi-\gamma\right)}$ against the
arrow. Under the basis $\{\mathbf{a}_{1},\mathbf{a}_{2}\}$ shown
in Fig. \ref{sfig:numerical-comparison-winding}(c), the Hamiltonian
reads, 
\begin{align}
H= & m\sum_{r_{1},r_{2}}\left(c_{r_{1},r_{2},\mathrm{a}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}-c_{r_{1},r_{2},\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{b}}\right)+\nonumber \\
 & +t_{1}\sum_{r_{1},r_{2}}\left(c_{r_{1},r_{2},\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+c_{r_{1},r_{2}+1,\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+c_{r_{1}-1,r_{2},\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+\text{h.c.}\right)+\nonumber \\
 & +t_{2}e^{i\gamma}\sum_{r_{1},r_{2}}\left[e^{i\varphi}\left(c_{r_{1}-1,r_{2},\mathrm{a}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+c_{r_{1},r_{2}-1,\mathrm{a}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}+c_{r_{1}+1,r_{2}+1,\mathrm{a}}^{\dagger}c_{r_{1},r_{2},\mathrm{a}}\right)+\right.\nonumber \\
 & \left.+e^{-i\varphi}\left(c_{r_{1}-1,r_{2},\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{b}}+c_{r_{1},r_{2}-1,\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{b}}+c_{r_{1}+1,r_{2}+1,\mathrm{b}}^{\dagger}c_{r_{1},r_{2},\mathrm{b}}\right)+\text{h.c.}\right],
\end{align}
where $c_{r_{1},r_{2},\mu},r_{1}\in\mathbb{Z},r_{2}\in\mathbb{Z},\mu=\mathrm{a},\mathrm{b}$
is the annihilator at the sublattice $\mu$ in the unit cell with
coordinate $r_{1}\mathbf{a}_{1}+r_{2}\mathbf{a}_{2}$, and $m\in\mathbb{R}$
is the detuning between the site $\mathrm{a}$ and the site $\mathrm{b}$.
In the numerical calculation, the parameters $t_{1},t_{2},\gamma,\varphi$
and $m$ are randomly generated as $t_{1}=0.70502$, $t_{2}=-1.32760$,
$\gamma=2.15618$, $\varphi=0.05877$, and $m=-0.64569$. As shown
in Fig. \ref{sfig:numerical-comparison-winding}(d), we calculate
$W(E,r)$ and $W_{\text{strip}}(E,r)$ against $r$ in the direction
of $\mathbf{a}_{1}$ for four random reference energies, which are
(i) $0.01162+0.68736i$, (ii) $0.54100-1.99811i$, (iii) $0.59046+1.89903i$
and (iv) $0.10591-0.34594i$. For the quasi-1D winding number $W(E,r)$,
$L_{2}$ is set to $20$ (that is, $40$ sites in the supercell).
Similarly to the case of the 2D HN model, the curves of $W(E,r)/L_{2}$
fit well with $W_{\text{strip}}(E,r)$ for the four samples.

\section{Amoeba spectrum of the 2D HN model}

\label{app:calculation-Amoeba}

In this section, we will prove that the Amoeba spectrum of the 2D
HN model is the same as the $x$($y$)-SGBZ spectrum, i.e., 
\begin{equation}
\sigma_{\text{Amoeba}}=\left\{ E\in\mathbb{C}\mid E=2e^{i\delta_{x}}\text{Re}\left(J_{x}^{*}e^{i\theta_{x}}\right)+2e^{i\delta_{y}}\text{Re}\left(J_{y}^{*}e^{i\theta_{y}}\right),\theta_{x},\theta_{y}\in\left[-\pi,\pi\right]\right\} ,\label{seq:Amoeba-spectrum}
\end{equation}
shown as Fig. \ref{sfig:calculation-Amoeba}(a). To prove Eq. (\ref{seq:Amoeba-spectrum}),
we need to prove that Amoeba $\mathcal{A}(E)$ exhibits a central
hole when $E$ is outside the spectrum and has no central holes when
$E$ is inside the spectrum.

When $|\beta_{x}|=\exp(\mu_{x})$ and $|\beta_{y}|=\exp(\mu_{y})$,
the ChP reads, 
\begin{align}
f_{xy}\left(E,e^{\mu_{x}+i\theta_{x}},e^{\mu_{y}+i\theta_{y}}\right) & =E-e^{i\delta_{x}}\left(\gamma_{x}e^{-\mu_{x}}J_{x}e^{-i\theta_{x}}+\gamma_{x}^{-1}e^{\mu_{x}}J_{x}^{*}e^{i\theta_{x}}\right)-e^{i\delta_{y}}\left(\gamma_{y}e^{-\mu_{y}}J_{y}e^{-i\theta_{y}}+\gamma_{y}^{-1}e^{\mu_{y}}J_{y}^{*}e^{i\theta_{y}}\right),\nonumber \\
 & =E-U_{x}\left(\theta_{x}\right)-U_{y}\left(\theta_{y}\right),
\end{align}
where the functions $U_{j}(\theta_{j}),j=x,y$ are defined as, 
\begin{align}
U_{j}\left(\theta_{j}\right) & \equiv e^{i\delta_{j}}\left(\gamma_{j}e^{-\mu_{j}}J_{j}e^{-i\theta_{j}}+\gamma_{j}^{-1}e^{\mu_{j}}J_{j}^{*}e^{i\theta_{j}}\right),\nonumber \\
 & =e^{i\delta_{j}}\left[\left(\gamma_{j}e^{-\mu_{j}}+\gamma_{j}^{-1}e^{\mu_{j}}\right)\text{Re}\left(J_{j}^{*}e^{i\theta_{j}}\right)+i\left(\gamma_{j}^{-1}e^{\mu_{j}}-\gamma_{j}e^{-\mu_{j}}\right)\text{Im}\left(J_{j}^{*}e^{i\theta_{j}}\right)\right].
\end{align}
The trajectory of $v_{j}(\theta_{j})$ is an ellipse on the complex
plane, whose major axis is parallel to $\exp(i\delta_{j})$.

By definition, $(\mu_{x},\mu_{y})$ belongs to $\mathcal{A}(E)$ when
and only when the ChP $f_{xy}(E,e^{\mu_{x}+i\theta_{x}},e^{\mu_{y}+i\theta_{y}})$
vanishes for some $\theta_{x}$ and $\theta_{y}$, or equivalently
when the ellipse $E-U_{x}(\theta_{x})$ intersects with the ellipse
$U_{y}(\theta_{y})$. Figure \ref{sfig:calculation-Amoeba}(b) and
\ref{sfig:calculation-Amoeba}(c) shows the cases for $E_{1}\notin\sigma_{\text{Amoeba}}$
and $E_{2}\in\sigma_{\text{Amoeba}}$, respectively. According to
the geometric-arithmetic mean inequality, the length of the major
semi-axis is $(\gamma_{j}e^{-\mu_{j}}+\gamma_{j}^{-1}e^{\mu_{j}})|J_{j}|\geq2|J_{j}|$.
Therefore, the ellipse $U_{j}(\theta_{j})$ must enclose (or coincide
with) the segment $2e^{i\delta_{j}}\text{Re}(J_{j}^{*}e^{i\theta_{j}}),\theta_{j}\in[-\pi,\pi]$,
shown as the dotted segments in Fig. \ref{sfig:calculation-Amoeba}(b)
and \ref{sfig:calculation-Amoeba}(c).

For $E_{1}\notin\sigma_{\text{Amoeba}}$, as shown in Fig. \ref{sfig:calculation-Amoeba}(b),
the two dotted segments have no intersections. When $(\mu_{x},\mu_{y})$
is close to $(\ln\gamma_{x},\ln\gamma_{y})$ enough, the ellipses
of $E_{1}-U_{x}(\theta_{x})$ and $U_{y}(\theta_{y})$ are separated
to each other, so that $(\mu_{x},\mu_{y})\notin\mathcal{A}(E_{1})$
in this case. For the derivative of the Ronkin functions, both the
winding number of $E_{1}-U_{x}(\theta_{x})$ around $U_{y}(\theta_{y})$,
and the winding number of $U_{y}(\theta_{y})$ around $E_{1}-U_{x}(\theta_{x})$,
vanish for arbitrary $\theta_{x}$ and $\theta_{y}$, so that the
derivative of the Ronkin function vanishes in this case. Therefore,
the Amoeba $\mathcal{A}(E_{1})$ shows a central hole in the neighborhood
of $(\ln\gamma_{x},\ln\gamma_{y})$.

For $E_{2}\in\sigma_{\text{Amoeba}}$, as shown in Fig. \ref{sfig:calculation-Amoeba}(c),
the two dotted segments intersect each other. Therefore, the two ellipses
of $E_{2}-U_{x}(\theta_{x})$ and $U_{y}(\theta_{y})$ do not have
intersections only when one ellipse encloses the other. When the two
ellipses intersect, $(\mu_{x},\mu_{y})\in\mathcal{A}(E_{2})$. When
one ellipse encloses the other, the winding number of the outer ellipse
around the points on the inner ellipse is non-zero, resulting in the
non-zero gradient of the Ronkin function. Consequently, $\mathcal{A}(E_{2})$
does not have central holes.

\begin{figure}
\centering

\includegraphics{Figures/sfig-calculation-Amoeba}

\caption{Sketches for the Amoeba spectrum of the 2D HN model. (a) Amoeba spectrum
of the 2D HN model. $E_{1}$ and $E_{2}$ are two reference energies
out of and in the spectrum, respectively. (b, c) Sketches of the ellipses
$E-U_{x}(\theta_{x})$ and $U_{y}(\theta_{y})$ for (b) $E=E_{1}$
and (c) $E=E_{2}$. The black dotted segments enclosed by the green
and orange ellipses denote the trajectories of $E-2e^{i\delta_{x}}\text{Re}(J_{x}^{*}e^{i\theta_{x}})$
and $2e^{i\delta_{y}}\text{Re}(J_{y}^{*}e^{i\theta_{y}})$, respectively.}
\label{sfig:calculation-Amoeba} 
\end{figure}

\bibliography{2D-skin-effect}
\end{document}

